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Non-thermal dark matter generation is an appealing alternative to the standard paradigm of 
thermal WIMP dark matter. We reconsider non-thermal production mechanisms in a systematic 
way, and develop a numerical code for accurate computations of the dark matter relic density. We 
discuss in particular scenarios with long-lived massive states decaying into dark matter particles, 
appearing naturally in several beyond the standard model theories, such as supergravity and su- 
perstring frameworks. Since non-thermal production favors dark matter candidates with large pair 
annihilation rates, we analyze the possible connection with the anomalies detected in the lepton 
_ cosmic-ray flux by Pamela and Fermi. Concentrating on supersymmetric models, we consider the 

>— ^ , effect of these non-standard cosmologies in selecting a preferred mass scale for the lightest super- 

symmetric particle as dark matter candidate, and the consequent impact on the interpretation of 
5-H ' new physics discovered or excluded at the LHC. Finally, we examine a rather predictive model, the 

G2-MSSM, investigating some of the standard assumptions usually implemented in the solution of 
■ the Boltzmann equation for the dark matter component, including coannihilations. We question 

the hypothesis that kinetic equilibrium holds along the whole phase of dark matter generation, and 
the validity of the factorization usually implemented to rewrite the system of coupled Boltzmann 
equation for each coannihilating species as a single equation for the sum of all the number densities. 
As a byproduct we develop here a formalism to compute the kinetic decoupling temperature in case 
(~ l \ of coannihilating particles, which can be applied also to other particle physics frameworks, and also 

^ |. to standard thermal relics within a standard cosmology. 
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I. INTRODUCTION 



The identification of the dark matter (DM) component of the Universe is one of the most pressing issues in Science 
today. Cosmological and astrophysical observations give compelling evidence for DM on a very wide range of scales; 
\l ' on the other hand, from a particle physics perspective, these observations have not provided any clear indications on 
the relevant properties of DM particles, such as their mass and the interaction strength with ordinary matter (for a 
recent review on the DM problem and its particle physics implications, see, e.g., (H). Viable DM candidates proposed 
in the literature include particles with a mass (close to) the Planck scale that are only gravitationally interacting, see 
e.g. 0, as well as ultra-light scalar particles possibly forming a condensate, see e.g. Q- The picture becomes more 
constrained once a mechanism to generate the DM term is considered. One of the most popular scenarios relies on a 
very general and elegant argument: stable massive species have an early- Universe thermal relic abundance scaling with 
d • the inverse of their pair annihilation rate, matching the cosmologically measured DM density when the annihilation 
cross section is about 3 • 10^^^ cm^s^^, a natural value for weak- force type couplings. This is the well-celebrated 
WIMP (weakly-interacting massive particle) miracle, allowing to embed a DM candidate in most of the proposed 
extensions to the standard model (SM) of particle physics, such as with the lightest neutralino in R-parity conserving 
supersymmetric (SUSY) extensions to the SM, or a heavy photon in T-parity conserving versions of Little Higgs 
models Q. One aspect which is particularly appealing is the fact that, in most of these examples, the existence of 
such states and of the symmetry enforcing their stability are not properties introduced ad-hoc to address the DM 
problem, but rather a byproduct of other features in the theory. 

Although extremely successful and attractive, the WIMP scenario faces also a few shortcomings. One the most 
severe is the fact that the idea of thermal generation of DM is based on the extrapolation of the properties of Universe 
from the earliest epoch at which the standard model for cosmology is well-tested, the onset of the synthesis of light 
elements (Big Bang Nucleosynthesis, BBN) at a temperature Tbbn of about 1 MeV, to the much earlier epoch of 
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WIMP thermal freeze-out, at the temperature Ttj.o. of about one twentieth of the WIMP mass. The WIMP miracle 
relies on three main assumptions: i) at freeze-out the Universe is in a radiation dominated phase with effective 
number of relativistic degrees of freedom as inferred from the SM particle spectrum; ii) there is no entropy injection 
intervening between Ttj.o. and Tbbn] Hi) there is no extra source of DM particles on top of the thermal component. 
There are particle physics models in which all these three hypothesis are actually strongly violated, such as in any 
theory containing heavy states that are very weakly (e.g., gravitationally) coupled to ordinary matter, such as the 
gravitino or moduli fields in SUSY setups. These states do not thermalize in the early Universe, they may dominate 
the Universe energy density, they are long-lived and potentially a copious source of entropy and DM particles at decay. 
The so-called cosmological gravitino or moduli problem refers to the very severe observational limits one encounters 
when these phenomena intervene during or after the BBN; on the other hand they can be perfectly consistent with 
available data if the lifetime of these fields is shorter than the age of the Universe at the onset of BBN, about 1 s, or, 
equivalently, if the Universe is reheated to a temperature Trh larger than Tbbn, where the reheating temperature 
Trh is defined as the temperature at which the Universe starts evolving in according to a radiation dominated phase 
after the field decays. 

The prediction for the relic density of DM in case of long-lived heavy fields is generally very model dependent; there 
are however a few definite pictures. One of the most appealing is, e.g., the one pointed in Ref. [5!] for SUSY DM in the 
anomaly-mediated SUSY breaking framework: it foresees the existence of a heavy modulus driving the Universe into 
a matter dominated phase and then decaying with a very large entropy production. The entropy injection reheats the 
Universe to a Trh in the range between a few MeV and about 100 MeV, and dilutes thermal relics (and gravitinos); in 
case the decay branching ratio into DM particles is unsuppressed, a very large number of DM particles is produced as 
well but gets instantaneously reduced to level at which their pair annihilation rate matches the value of the Universe 
expansion rate at Trh- As a rule of thumb, the relic density of DM in this scenario turns out to be approximately equal 
to the result for thermal relic WIMPs scaled up by the ratio Ttj.o./Tjm: to be in agreement with the cosmological 
measurements of the DM density, the DM pair annihilation cross section needs also to increase accordingly (this 
solution to the moduli problem and the new match in the annihilation cross section appear naturally in some particle 
physics frameworks and is sometimes referred to as a 'non-thermal WIMP miracle' [6|; a more detailed discussion of 
the framework and a list of relevant references is given below). 

The interest in DM particles with large pair annihilation cross section has been recently triggered by the fact 
that cross sections larger than the standard face-value for WIMP DM are needed to provide a DM positron source 
accounting for the rise in the positron fraction in the local cosmic-rays measured by the PAMELA detector 
The picture of non-thermal generation of DM has however a much broader phenomenological impact, e.g., shifting 
significantly the mass scale for which a DM particle embedded in a SM extension is cosmologically relevant or excluded, 
with a direct impact as well on the interpretation of new physics discovered or excluded at the LHC. 

In the paper we reexamine the issue of non-thermal generation of DM implementing a full numerical solution of the 
relevant set of equations, including the equation of motion for the heavy fields and the system of coupled Boltzmann 
equation, and avoiding to introduce approximations such as instantaneous reheating and production of DM particles, 
or others. The code is interfaced to an appropriately modified version of the public available DarkSUSY package 
and allows, for any definite particle physics scenario, a very accurate computation of the relic abundance of the 
DM particle. Examples of its applications are given in a sample of progressively refined particle physics scenarios: 
we first consider a toy-model in which DM particles are schematically defined through the values of the mass and 
a temperature-independent pair annihilation cross section. We examine then the impact of heavy long-lived fields 
on the phenomenology of the most widely studied case for WIMP DM, i.e. the case of neutralino dark matter in 
the minimal supersymmetric extension of the standard model (MSSM), focussing in particular on the framework 
usually dubbed " Split Supersymmetry" [lo| , in which, since the sfermion sector is not playing any relevant role, 
the parameter space is sharply reduced compared to a g eneral MSSM, but it is still general enough for our purpose. 
Finally we examine in detail the G2-MSSM scenario , a particular class of theories with rather precise predictions 
for the spectrum of low energy SUSY particles as well as of the gravitino and the moduli fields; in such a definite 
framework it is interesting to test the validity of some of the approximations that are usually given for granted in the 
solution of the Boltzmann equation for DM: In this model the lightest neutralino is the lightest SUSY particle but it 
also nearly degenerate in mass with a chargino; coannihilations, namely the interplay between the two in the early 
Universe, are usually treated by writing a single Boltzmann equation for the sum of the two species, however for very 
low reheating temperatures this approach may not be valid. We solve the full system of coupled Boltzmann equation 
and address also the issue of energy losses for relativistic neutralinos and charginos injected by moduli decays. As 
a byproduct and final step, we develop here for the first time the formalism to compute the DM kinetic decoupling 
temperature for a system undergoing a low-temperature reheating and for which kinetic equilibrium is maintained in 
a chain of coupled processes rather than by the elastic scattering of a single DM particle on thermal bath particles. 
This is an important result since thermal kinetic decoupling, namely when DM scattering goes out of equilibrium 
(as opposed to the thermal chemical decoupling mentioned above which refer to the departure from equilibrium of 
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the pair annihilation processes) determines the small-scale cutoff in the spectrum of matter density fluctuations, see, 

The paper is organized as follows: In Section 2 we review the particle physics scenario and illustrate the approach 
we follow to trace the evolution of the long-lived heavy fields and compute the relic density of DM particles. In 
Section 3 we give a first example of this procedure introducing a simple toy-model; for this model we cross check 
various approximate scalings for the DM relic density discussed in the literature, and also derive, under different 
configurations, what range of reheating temperature would be needed to provide an explanation for the very large 
annihilation cross section which may be inferred from the PAMELA positron excess. In Section 4, still within a 
schematic treatment of the decaying heavy fields, we consider the case of split SUSY in the MSSM and discuss the 
impact of low reheating temperatures on the cosmologically relevant portion of its parameter space. In Section 5, we 
focus on the G2-MSSM, with the numerical treatment of the sequence of steps induced by the set of heavy long-lived 
fields present in the model; the treatment of this case is refined in Section 6, where we also compute the DM kinetic 
decoupling temperature for this model. 



II. THE GENERAL FRAMEWORK FOR NON-THERMAL DARK MATTER PRODUCTION 

We consider a particle physics framework embedding three sets of beyond-SM states: Let {xa} be a set of particles 
sharing a conserved quantum number and having sizable couplings with SM particles; the first property ensures that 
the lightest of them, say xoi is stable, while the second guarantees thermalization at sufficiently large temperatures. 
We also require that xo has zero electric and color charges, so that it can play the role of DM candidate. The 
fields {Xi} in the second set, which we will refer to as cosmological moduli, are instead heavy states which are very 
weakly interacting, out of thermal equilibrium in the early Universe and long-lived; we assume they can condensate, 
potentially dominate the Universe energy density at an intermediate stage in its evolution, and later decay producing 
both SM particles, with a sharp increase in the entropy density, and Xa fields. Finally the framework may contain also 
additional long-lived states, say ipi, out of thermal equilibrium but with a subdominant contribution to the energy 
density, possibly sharing the quantum number protecting the stability of xo ! these particles may also be produced in 
the decay of the Xi states. 

A typical scenario of this kind is provided by SUSY extensions of the SM. We will consider in particular cases 
in which the states {xa} are the superpartners of SM particles within the MSSM, and R-parity is the symmetry 
protecting the stability of the lightest SUSY particle (LSP), which we will assume to be the lightest neutralino. For 
what regards the moduli {^i}, there are several possibilities. First of all, it is quite common, in SUSY theories, to 
find field configurations for which the scalar potential is fiat; these configurations are referred as 'flat directions' and 
can be described by a chiral superfield. For our purposes, only the scalar component of these multiplets is relevant; we 
will refer to this part with the term modulus. SUSY breaking can lift the flat directions inducing a mass term for the 
moduli in the scalar potential. Further candidates for long-lived states arise in supergravity theories. The gravitino 
is not in thermal equilibrium in the early Universe, however, plays a different role compared to moduli fields; we will 
restrict to the case when gravitinos are heavy and not the LSP (otherwise the phenomenology would be very different 
from the one discussed in this paper), falling in the category of the i/ji fields introduced above. Another possibility is 
the Polonyi field [H-Ell which is introduced in many SUSY breaking schemes. Finally supergravity can be seen as a 
low energy limit of string theory, in which scalar fields can appear in the compactification of extra dimensions. 

The impact on cosmology of the moduli can easily sketched under a few simplifying assumptions. Consider for 
simplicity the cases of a single modulus X which decays when it dominates the Universe energy density; focussing on 
Planck suppressed interactions, its decay width can be written in the form: 

^-^^-mI' W 

with Dx some coefficient depending on the specific model, mx the field mass and Mpi the reduced Planck mass. 
Assuming instantaneous conversion of the energy density into radiation, one usually defines the reheating temperature 
Trh through the expression: 



= 90 M^i' 



where gce(T^ii) is the effective number of relativistic degrees of freedom at Trh- Inverting this expression, one finds 
approximately that the onset of the standard radiation dominated phase happens at the temperature: 



IMeV ■ ^ [geff(T] 



^RH 
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To avoid spoiling predictions of the standard BBN, one needs to require Trh ^ 4 MeV [it'], which, for Dx of order 
one, translates into a lower limit on the mass of the cosmological modulus of about 30 TeV. 

At this level of approximation the evolution of the system would be fully specified by the X decay width (or Trh) 
and the amount of energy density converted into dark matter particles (which, in the treatment above, was implicitly 
assumed to be tiny compared to amount going into radiation). Having instead in mind to be able to treat a system 
in which, for the full set of Xi fields, spectra, lifetimes and branching ratios in the decay are calculable in the given 
particle scenario, we will not refer to instantaneous reheating, but rather follow more explicitly the evolution of the 
moduli. In principle this could be done by studying a full set of coupled equations of motion, having specified the 
potentials for each field [HI, [l^. The result would be that, at early times, each field stays frozen in a time dependent 
minimum; when becomes of the order of the Xi mass, the equation of motion takes the form of the one for a 
damped harmonic oscillator. The oscillations satisfy a pressureless equation of state and hence the scalar field behaves 
like a condensate evolving as a matter fluid; provided that enough energy is initially stored in Xi , the Universe enters 
a phase of matter domination lasting until the field decays, with the transition that needs to be treated as a continuos 
process. In practice, even for models for which the physics related to the Xi fields is given in some detail, it is 
difficult to describe potentials and their temperature evolutions beyond the toy model level; for our purposes it will 
be sufficient to follow the evolution of the system starting from the phase of coherent oscillations. Each state Xi is 
then traced through an equation for its energy density: 

^ + 3Hpx, = -Tx^px, , (4) 

and in case of several moduli present at the same time, the single equations are included in the system at the time 
t — 3/2H — mxi, assuming the energy density stored in the field at this time is equal to (1/2)771^ Mp,[T8l-l20j: we 
will comment later on the fact that the final density of dark matter is not sensitive to these assumptions. 

The decay of the Xi particles produces SM particles, Xa states cascading to the DM particle, and, eventually, 
the long-lived ijji fields, in turn decaying into radiation and, possibly, DM particles. From the first principle of 
thermodynamics, one can write an equation for the total energy density and pressure associated to SM, Xa and tpj 
states, respectively, p and p, in an implicit form: 

^ + 3i/(p + p) = ^rx.px.. (5) 

i 

This equation is treatable once separating p and p in components. Starting with the ipj particles, one can safely assume 
that they are produced in given number at decays, get diluted and redshifted by the Universe expansion without 
interacting with other species and decay themselves (an eventual term associated to the production via inelastic 
scattering off SM or Xa particles is not introduced since such term becomes relevant only at large temperatures, while 
we will only consider here the case of moderate to low Trh; also, we are not considering the possibility of a "0^ particle 
decaying into a lighter i/jfe state, since we will not encounter a case of this kind in explicit models and it would just 
complicate the notation). The Boltzmann equation for the ijjj number density is: 

— + iHn^^ = > — i Tx^px, - r^ .n^^ , (6) 

at ^-^ mx 

i * 

where B^ .^Xi is the mean number of particles ipj generated in the decay of the field Xi, i.e. it is the product of the 
branching ratio of decay into ipj times the mean multiplicity. 

To trace the number density of the Xa states, especially when two or more of these are nearly degenerate in mass 
(coannihilating particles), one should refer to a system of coupled Boltzmann equations describing: the source from 
the decay of the Xi and ipj fields; their changes in number density due to pair production from and annihilation 
into SM particles; the energy exchanges with SM thermal bath particles through elastic scatterings processes; the 
redistribution in the relative number density by inelastic scattering of a given Xa into a different Xb state; the decays 
of Xc into lighter Xd particles and and of these to the lightest stable species. This is usually not done since it is a 
system of coupled stiff equations one needs to solve numerically; moreover it is usually not necessary to do it, since 
one is interested only in the number density of the lightest state after all heavy states have decayed into the stable 
one. Rather than tracing the number density n-^^ of the individual state Xa, one usually solves a single equation 
written for the sum of all the number densities, — X^a^-Xa; i-^- 

^+3Hn^^ -(aeff^.) [n^ - (r^^^f] + ^ ^Tx.Px, + ^ i?^T^,n^, (7) 
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where = J2a ^Xa,Xi and B.^^ = J2a^Xa,tl'j have been defined in analogy to B^.^Xi- In this equation stands 
for the sum of thermal equilibrium number densities, and the term proportional its square accounts for the production 
of particles Xa in pair annihilations of SM thermal bath particles, while the effective thermally averaged annihilation 
cross section: 



a.b 



Xb 
eq 



(8) 



is written as a weighted sum over the thermally averaged annihilation cross section of any Xa-Xb pair into SM particles; 
the processes giving a sizable contribution to this sum are only those for which the mass splitting between a state Xa 
and the lightest state xo are comparable to the thermal bath temperature T . There are two main assumptions which 
allow to implement Eq. ([T]) to trace n^. The first is kinetic equilibrium for each species Xaj namely that the scattering 
processes on thermal bath particles are efficient and make the phase space densities for each particle trace the spectral 
shape of the corresponding thermal equilibrium phase space density, namely fa{ka,t) — C{t) ■ f^'^{ka,t) (with the 
coefficient C depending on time but not on momentum) . Within this assumption, we treat as instantaneous the energy 
depletion from the relativistic regime when particles are injected from moduli or ipj decays to the non-relativistic 
velocities in the low reheating temperature background plasma; it also allows to factorize the number density out 
of each thermally averaged annihilation cross section (which is defined in terms of thermal equilibrium phase space 
densities). To implement the factorization of the individual terms in the sum of Eq. ([8|) one needs also to assume 
that n^ci/uy^ ~ n'^^/n'!^, a quantity which, in the Maxwell-Boltzmann approximation for the equilibrium phase space 
densities, as appropriate for non-relativistic particles, is proportional to the number of internal degrees of freedom 
g^^ and is exponentially suppressed with the ratio between mass splitting and temperature; this approximation is 
strictly valid only in case inelastic scatterings of Xa particles are efficient over the whole time interval in which the 
pair annihilation term is relevant. Within the standard computation of the thermal relic density for WIMPs, the two 
assumption are in general well justified, since the kinetic decoupling and the decoupling of inelastic scatterings usually 
take place at a much lower temperature than chemical decoupling; while assuming that Eq. ([T]) is valid in the next 
Sections, in Section 6 we study this issue in more details and, considering a specific particle physics scenario, address 
the problem of kinetic decoupling in models with non-thermal generation of DM particles. 

We keep track of the SM states only through their contribution to the radiation energy density and pressure which, 
using Eq. ([5]) and subtracting the contribution from ipj and Xa fields, obey the equation: 

— +ZH{pr+pr) ~ 2^11 \Vx,px,+2_^[{E^,) ~m^B^^)T^^n^. 

-fm^(aefft.)[n2-(n^^)2] . (9) 
In this equation {E^.^Xi) is the mean energy of the particle V'j at injection from the decay of the modulus Xi: 



dE' 



dE' 



(10) 



with dhf^-.Xi/dE' the energy spectrum from the decay normalized to 1; (E^.) is instead the mean energy for ipj 
particles: 



fdt'j:^Tx.PxAt') fdE''^^-''^ 
Jo , mx, J 



dE' 



-,1/2 



(11) 



■exp[-rv,,(t-i')] 



Finally, in Eq. ^ we have assumed that the mean energy of the Xa states is equal to the mass of the lightest 
state xoy neglecting, at this level, thermal corrections and mass splittings between the coannihilating states, as well 
as the pressure term associated to Xa- 

Eqs. dH), dH), and ^ define a system of coupled equations, closed by Friedmann equation giving H. In its 
numerical solution, it is more convenient to use as independent variable, rather than the time t, the rescaled scale 
factor A = a/aj, with aj an arbitrary parameter with dimension of the inverse of an energy. Following [20j . we will 
use as dependent variables the dimensionless quantities: 



^x. 



PXja 
A 



N,!, 



and 



Ny EE n^a^ 



(12) 
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with A an arbitrary energy scale, plus the temperature T, expressing pji and pn in terms of the entropy density 
through the standard definitions: 



s{T) 



Pr{T)+pr{T) _ 2tt 



T 



45 



^/.eff(T)T3 and p„(T) ^^g,^{T)T' = --^p^^Ts{T) ^ (13) 



3 5cff(r), 



4/leff(T)^ 



with QcS and hcS the effective number of relativistic degrees of freedom. The values of a/ and A are chosen in order 
to guarantee the best numerical stability to the solution, a sample guess being, respectively, and Trh, with the 
approximate reheating scale as given through Eq. ([2]). After this change of variables the system becomes: 



dA 



dA 



A^I-'aT 

A^l^aT 
H 



Xi 



B 



(14) 



dN^ 
dA 

dT 
dA 



A^/^a/^n 



A^r^aT 



= 1 



T dg^ff 



4geff dT 



hcs T 



^eff 



5off A 35cffs(r) ^5/2^3/2^ 



.Xi 



mxi 



^x,i 



Xi 



where H is defined from the Universe expansion rate, as: 



n = {aiAf'H = 



1/2 



(15) 



The relic density of dark matter can be evaluated by evolving these equations from an initial time, which we assume 
to be the time when the heaviest modulus starts its coherent oscillations, up to the stage when the DM comoving 
number density becomes constant. 



III. NON-THERMAL DM PRODUCTION IN A TOY MODEL AND RELEVANCE FOR PAMELA 

We discuss first a minimal framework with a single cosmological modulus X decaying into the DM particle x- 
Rather than detailing a specific particle physics scenario, in this first example we define x only through its mass 
and pair annihilation rate into SM particles, whose thermal average is assumed not to depend on temperature, as 
appropriate for S-wave annihilations. We also avoid dealing with eventual other states charged under the quantum 
number protecting the stability of Xi assuming that they have a sizable mass splitting with respect to x, and hence 
have very short lifetimes and do not enter in the Boltzmann equation for x- Under these hypotheses the system of 
coupled equations reduces to three equations only: the first for the decaying modulus, the second for number density 
of particle Xi sourced from the decay and depleted by pair annihilations, and the last for the temperature. 

In this simplified picture, the main trends in the non-thermal DM production can be illustrated even at the level 
of approximate analytical formulae; we briefly summarize here some of these features, as we will recover them in the 
numerical solution of this model as well as in the more involved scenarios we will consider later (for a more detailed 
discussion, see, e.g., [111). First of all, if the modulus decay induces a large increase in the entropy density and this 
happens at a later stage with respect to the chemical decoupling for x, the thermal relic density of x is greatly diluted 
and can be neglected, with the only relevant x source being the particles produced in the decay itself. The entropy 
injection is a continuos process making the reheating phase last for an extended period during which one can show that 
the temperature evolves as T cx a~^/^ and the universe expansion rate as iJ oc T"* (20l.[23|. A standard approximation 
is however to treat the decay of the field and the thermalization of the products as instantaneous processes, and define 
the reheating temperature Trh according to Eq. ([2]); depending on whether at Trh the dark matter pair annihilation 
rate T = n^{<Jv) is larger or smaller than the expansion rate H, there are two distinct regimes determining the relic 
density for x |5l, l24j . If F is much larger than H, pair annihilations are very efficient and instantaneously decrease in 
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the number density of % to the critical density level corresponding to F ~ when the annihilations stop; such critical 
density is then simply equal to: 

^ (16) 



iav) 



As usually done, we normalize the number density to the entropy density introducing the quantity Y-^ = '^xl ^i since 
when annihilations become inefficient, if there are no further entropy injection phases, such ratio becomes constant 
and can be used to estimate the relic density for 



^^X - 77^''x\^0) - TT^T^^xl-'RHj OC — , (17) 



where Pc{Tq) and s(To) refer to the Universe critical density and entropy density at present. The rule of thumb F ~ 
is the same criterium implemented for an approximate estimate of the relic density in the standard thermal decoupling 
picture for WIMP dark matter, except that the reference temperature in this latter case is the thermal freeze-out 
temperature Tt.f.o. — m^/20. Following the same steps, one finds that the thermal relic density fij scales with the 
inverse of Ttj.o., and hence that the relations of i^^^ with 57j and the WIMP pair annihilation cross section are 
approximately given by: 



,NTu2 ^ ^T^2 ^ ^x/20 _ 3 • 10-26 cm^ s-1 

Trh ^ Trh 



A particle x whose thermal relic density is small compared to the DM density because the annihilation rate is too 
large, may become a viable dark matter candidate for an appropriate value of Trh- This simple rescaling holds 
whenever the particles x ^re copiously produced in the modulus decay and if the pair annihilation rate is sufficiently 
large; in the following, we refer this scenario as 'reannihilation regime'. If instead ri^(TRH) is lower than n^(TRH), 
the particles produced in the decay do not interact further and their number density per comoving volume is frozen, 
being: 



s(Trh) mx s(rRH) 4tox 



and hence giving a non-thermal relic density which is about (see also, e.g., |23|): 



n^^h^ ^ 0.2 . W^Bx^^^ (20) 

^ mx IMcVlOOGcV ^ ' 



Note that, in this case, the final dark matter density depends on the physics of moduli not only through its propor- 
tionality to the reheating temperature but also through the ratio between the average number of particles x produced 
per decay and the modulus mass Bx / mx ■ This non-thermal scaling applies to the cases in which either the pair 
annihilation rate is small or the average number of particles x produced per decay Bx is small. ^ 

We are now ready to discuss numerical results within this simplified scenario. As just outlined, the relevant 
parameters for the relic density calculation are the particle mass and pair annihilation cross section, as well as those 
setting the efficiency in producing dark matter particles and the energy density of the field at decay; regarding the 
latter we will treat as free parameters Bx and the mass of the modulus mx , which in turns sets the decay width Tx 
and hence the reheating energy (we start with the assumption of gravitational interactions in the decay, and comment 
shortly on how to interpret results in case of a more general expressions for Fx)- Since we are tracing the full evolution 
of the field X, we are not in the limit of instantaneous reheating and do not implement the definition of reheating 
temperature as quoted in Eq. the Trh we refer to when illustrating results is extrapolated from the numerical 
solution, matching the T cx a"^/^ scaling obtained in the phase when the X decays act as a large source of entropy to 
the T ~ scaling in the subsequent radiation dominated regime (this prescription of matching asymptotic solutions 
is not totally rigorous since we should also take into account eventual changes in the number of relativistic degrees 
of freedom contributing to the entropy density; in practice, however, since the transition between the two regimes is 



^ Ref. (23)1 classifies two extra scenarios, already studied, e.g., in Ref. [20il . corresponding to the case in which the main source of x particles 
is pair production from SM background states; these applies essentially only in the limit of Bx — > which we are not going to discuss, 
although the method outlined here would be suitable for them as well. 
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FIG. 1; The number density of the DM particles normalized the entropy density as a function of temperature, for few values 
of the reheating temperature and Bx fixed to 0.25 (left panel), and for a sample Trh while varying Bx (right panel). In the 
left panel we also plot with dashed-dotted lines the quasi-equilibrium density K^^^(r) for each Trh, and with a dashed line 



the critical density Y-^{T) for the case Trh = 
panel the dashed line shows Y-^{T) for Trh = 
section being, respectively, 1 TeV and 5 • 10" 



100 GeV in which reheating takes place before thermal freeze out. In the right 
180 MeV. Results refer to a DM particle with mass and pair annihilation cross 



always rather sharp, the Trh found in this way is always very accurate in parametrizing the total entropy injection 
from the X decay; note also that Trh is not used in any step of the numerical computation). 

In Fig. [T] we consider a sample DM particle with heavy mass, — 1 TeV, and large pair annihilation cross section, 
(av) — 5 • lO^^^'cm'^s"^; the ratio of the DM number density to the entropy density is plotted as a function of T 
(note that, since we want to compare directly Trh with T, rather than showing versus the inverse of temperature 
as usually done, we plot it versus T and use a logarithmic scale which decreases from left to right). In the left 
panel we have fixed Bx to a sample value representative of the case when the branching ratio of the decay into x is 
unsuppressed, and vary mx to select a few values of the reheating temperature; in the right panel, vice versa, we fix 
mx and vary Bx ■ The system of equations is solved assuming the initial energy density in the modulus is equal to 
1/2 rri^Mpi and that the radiation energy density is at the same level When Trh is larger than the thermal 

freeze-out temperature Ttj.o. for this model (the case for Trh — 100 GeV in the plot), the temperature evolution 
of Y^ is obviously the same as in a standard thermal WIMP framework: Y^ follows first the thermal equilibrium 
distribution along its Maxwell-Boltzmann tail, in a phase in which the main source of DM particles is pair production 
by SM background particles and this is balanced by DM pair annihilations, and then at Tt.f.o., when n'^ becomes 
smaller than and pair annihilations become inefficient, Y^ settles on a constant value. When Trh is reduced two 
effects intervene: first of all, the thermal freeze out temperature tends to increase since the modulus contribution 
to the Universe energy density increases H and hence n^; at the same time, the dominant source of DM particles 
becomes the modulus decays rather than SM pair creation. If x number density from the decay exceeds n^, this source 
term is balanced by DM pair annihilations and tracks the quasi-static equilibrium (QSE) density, as defined, e.g., 
in Ref M : 



^QSE ^ f Bx^XPX ^ 

^ \mx{(TcfLv) 



(21) 



For our sample DM model, this is the behavior we find in all cases with large Bx and Trh < Tt.f.o.'- starting at 
high T, Y^ follows first Y^'^ , then it becomes equal to Y^^^ up to about Trh when the modulus DM source drops 
exponentially, Y^^^ crosses and hence Y^ gets frozen. Regarding the temperature scalings in the plot, in the phase 
when the modulus dominates the energy density and is the main entropy source, we see that both Y^^^ and Y^ are 
proportional to T, except for a short low temperature phase in the examples for Trh = 30 and 6 MeV during which 
the entropy injected but the modulus decay is still negligible compared to the initial entropy and hence a oc T^^, 
making Y^^^ and Y^ rise as T^^l^ . For small Bx^ Y^^^ becomes smaller than 1?^, DM annihilations are inefficient 

X. K. A. A. 

and simply scales as Bx^xPx /"mx ■ t, up to the reheating temperature when the modulus source drops and Y-^ 
becomes constant; for what concerns the behavior in temperature, once again, in the phase in which the decay injects 
DM particles, the scaling just given translates into Y^ cx T, while for very small Bx one can also see a transient in 
which the amount of DM produced in the decay is small compared to the thermal component and Y^- simply refiects 
the entropy increase, decreasing faster than T. 
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FIG. 2: Left panel: scaling of the relic density with the parameter Bx for the same sample DM model selected for Fig [T] and 
for a few values of Trh- Right panel: relic density versus Trh, having fixed Bx = 10~^ and rescaling the value of (av). The 
black horizontal lines represent the cosmological DM density as extrapolated by the WMAP 7-year data [l^ ■ 



In the example displayed, the specific set of initial conditions implemented to solve the system of equations has a 
negligible impact on the final comoving density of DM particles. In fact the latter is insensitive to the choice of the 
initial energy density in the moduli and the relative weight with respect to the initial radiation energy density provided 
that the physical mechanism determining the DM relic density starts becoming efficient at temperatures lower than 
the temperature at which the scaling T ^ a~^^^ begins. More precisely: in all cases considered in this paper, the 
DM pair annihilation rate is large enough to guarantee, even in the non-standard cosmological scenarios considered 
here, chemical equilibrium at T > m^; the final relic densities is then determined by the physics taking place between 
the thermal freeze-out temperature and the reheating temperature. If the T ^ a^"^/^ scaling starts sufficiently earlier 
than Trh, the entropy production guarantees the suppression of the DM thermal component and, at the same time, 
variations in the entropy release with the field energy density are compensated by a different efficiency in the non 
thermal production, leaving then the final result unchanged. If, on the contrary, the DM thermal relic component 
is not totally diluted, the T ~ a~^^^ phase needs to start before the thermal freeze-out temperature, otherwise the 
variation of dilution due to entropy release stemming from the initial conditions has a direct impact on the relic 
density as well. 

In Fig. [5] we plot the relic densities for the x state. In the left panel we refer to the same model introduced for 
Fig. [TJ select a few values for the reheating temperature and display results as a function of Bx] as expected from 
the discussion above, one can see that ily, becomes essentially independent of Bx in the limit of large Bx, while it 
scales linearly with Bx when the modulus source function is too small to make exceeds n^. Also visible at large 
Bx is the scaling of the relic density with the inverse of Trh, as expected from the analytical estimate in Eq. (jl8l) . 
When Bx is small il^ is expected to scale with Trh/^x- In our approach Trh and nix are correlated; from the 
instantaneous approximation, Eq. one expects mx oc T^Q giving cx Tj^h^, which is approximately the scaling 
seen in the plot for very small Bx- The dependence on the reheating temperature is shown more explicitly in the 
right panel of Fig. [51 where, having fixed Bx to an intermediate value, we let the annihilation cross section vary of 
a few orders orders of magnitude around the value chosen for the plot on the left hand side; the relic density scales 
with the inverse of (av) whenever reannihilation takes place, while evidently the solution does not depend on (av) in 
case annihilation processes are inefficient. 

In most scenarios containing cosmological moduli it is hard to tune the model in such a way that very tiny Bx 
are obtained, hence the framework we are discussing becomes interesting mainly when x is associated to a large 
annihilation cross section, preventing the overproduction of DM with respect to the experimental bound. DM models 
with a (av) which is two or three orders of magnitudes larger than in the standard thermal relic scenario would 
be very interesting also from the point of view of indirect DM detection and have been invoked to address the 
excess in the lepton cosmic ray flux by Pamela and Fermi. In Fig. [3l choosing a few sample values of Trh and two 
representative cases for Bx, we scan the parameter space {{crv) - m^) searching for configurations in which the x 
relic density matches the central value for the cosmological DM density as estimated from the WMAP 7-year data, 
namely Vt^h? = 0.1123 ± 0.0035 A curve corresponding to a given Trh becomes horizontal when Trh becomes 
larger than Tg.f.o., i.e. we recover the standard thermal result of the relic density being independent of mass for 
S-wave annihilations; on the other hand it becomes vertical when annihilations become inefficient and hence 51^ stops 
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FIG. 3: Values of mx and (ctw) for which the rehc density of the % particles matches the cosmological DM density, for a few 
values of Trh and two representative cases for Bx- Also shown are the region in this plane compatible with the Pamela and 
Fermi electron/positron excesses in case of leptophilic DM 28], and the WMAP limit and Planck projected sensitivity [29l. |30|] 
stemming from the impact of residual DM annihilations on reionization. 



depending on (crw). In general going from large to small Bx-, keeping Trh fixed, shifts the results to smaller (ctw) 
and larger m^^. In the same plot, supposing we are now referring to a leptophilic DM candidate, namely annihilating 
democratically into the three lepton species (28j , we have superimposed the region in the parameter space which have 
been found to be compatible with the Pamela and Fermi electron and positron data, as derived, e.g., in Ref. psj : 
the comparison is meant to be qualitative since we are not considering here a detailed particle physics scenario, it 
shows however what are the main trends that should be fullfilled to find an agreement. Also shown is the bounds on 
leptophilic models following from WMAP CMB data [l^ [sOl : the limit stems from the impact of residual (namely 
much later than thermal decoupling) pair annihilations on reionization, and will be soon improved by the Planck 
experiment in case of no signal. 

The last issue we wish to discuss in this Section is an implicit dependence on the modulus mass mx we have 
ignored so far: As mentioned above we have been varying mx to retrieve different values of Trh as extrapolated from 
the numerical solution of the system of coupled equations; the underlying assumption here is that we computed the 
modulus decay assuming gravitational coupling and a two body final state. Having in mind more general scenarios 
like those, e.g., in Ref. [a. l3ll - [35| . we may consider replacing: 

„ 1 Tr^v m\- 

rx = ^^ ^ r^-A^' (22) 

where now Ae// encodes both the coupling of the effective operator responsible for the decay and the kinematical 
factors. From the approximation of instantaneous reheating one sees that, to keep Trh fixed after this replacement, 
one needs simply to approximately shift: 



A^ 



1/3 



rr^x ^ rnx-\-^] . (23) 



The modulus mass however appears explicitly also in Eq. ([7]) when, in the DM source function from modulus decays, 
one converts from the modulus energy density to its number density. To compensate for this and use results displayed 
in this and the next Sections, one then should also shift the values reference values for Bx as: 



2 \ 1/3 



IV. IMPACT OF A TOY MODEL ON THE MSSM 



In this Section, while still referring to the schematic picture with a single cosmological modulus X parametrized 
through its decay width Tx and the DM yield Bx, we introduce an explicit particle physics scenario for the {xa} fields. 
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considering neutralino DM in the MSSM. As already mentioned, a scenario with DM production from cosmological 
moduh can arise quite naturahy in supergravity/superstrings theories. Some scenarios, such as gauge mediated 
supersymmetry breaking are actually troublesome since the moduli tend to be light and decay after the onset of BBN, 
see, e.g. [1, H^-fssI] ; to solve this proljlem, one needs to invoke a mechanism of dilution of the moduli number density, 
making these models not viable for non-thermal DM production. 

In the MSSM there are four neutralinos, spin 1/2 Majorana fermions obtained as the superposition of two neutral 
gauginos, the Bino B and the Wino W^, and the two neutral Higgsinos and i?2 ■ 

X,^NaB + N,2W + NaH^ + N^iH^ , (25) 

where the coefficients Nij are the elements of the matrix which diagonalizes the neutralino mass matrix, and are 
mainly a function of the Bino and the Wino mass parameters Mi and M2, and of the Higgs superfield parameter fi, 
while depend rather weakly on tan /?, the ratio of the vacuum expectation values of the two neutral components of the 
SU(2) Higgs doublets, which appears in the off-diagonal terms of neutralino mass matrix. The hierarchy between Mi, 
M2 and fi sets hence whether the lightest neutralino, which we are assuming also as lightest supersymmetric particle 
(LSP) and stable, is mostly bino-, wino- or Higgsino-like. From the point of view of DM production in the early 
Universe, pure Binos have pair annihilation cross sections dominated by SM fermion final states, which are helicity 
suppressed (S-wave annihilation cross sections proportional to the square of the mass of the final state fermion) and 
scale approximately with the inverse of the forth power of the mass of the corresponding sfermions; given current 
accelerator bounds on sfermion masses, Binos tend in general to have a too large annihilation cross section to be 
thermal DM relics. On the other hand Winos and Higgsinos have unsuppressed annihilation cross sections into W 
bosons and their annihilation cross section tends to be too large for thermal production unless one considers heavy 
states, about 2.4 and 1.1 TeV, respectively, for a pure Wino and a pure Higgsino, since, in this case, the annihilation 
cross sections scales approximately with the inverse of the square of the neutralino mass (we will show results obtained 
computing tree-level annihilation amplitudes; in case of TeV Winos and Higgsinos actually the result would change 
slightly, shifting the masses to slightly larger values, when taking into account that, for such heavy states, the weak 
interaction acts as a long-range attractive force which deforms the wave function of the annihilating DM pair, an 
effect usually referred as Sommerfeld enhancement, see, e.g., and references therein). Since we will be mainly 
interested in discussing the shift on the mass scale for neutralino DM due to non-thermal effects, we refer here to a 
supersymmetric framework maximizing this dichotomy between underproduced and overproduced thermal states, the 
so-called " Split Supersymmetry" scenario 0, [13] • This indicates a generic supersymmetric extension to the SM in 
which fermionic superpartners have a low mass spectrum (say at the TeV scale or lower), while scalar superpartners 
are heavy, with a mass scale which can in principle range from hundreds of TeV up to the GUT or the Planck scale 
0, a feature which can occur in a wide class of theories, see, e.g. [iol - li^ . Leaving out of the discussion also gravitinos 
which are assumed to be heavy and not produced in the modulus decay, the system {xa} reduces to neutralinos and 
charginos, whose annihilation and coannihilation effects we treat interfacing the model the DarkSUSY package Q. 
Finally for what regards the Higgs sector, the scenario has only one light state SM-like Higgs; the value of its mass, as 
well as tan /3 have no sizable impact on the overall picture, hence we keep them fixed to sample values, respectively, 
114.4 GeV and 10. 

In Fig. m we scan the parameter space Mi, M2 and fi searching for models whose relic abundance matches the central 
value from the 7-year WMAP estimate of the DM density in the Universe. There are three pairs of plot in which we 
vary two of the parameters, fixing the third to a heavy scale; in each pair, one plot is for a large _Bx, while the other 
is for a small but not negligible Bx- As in the previous Section, we vary mx to change the reheating temperature 
scale, assuming a two-body gravitational decay for the modulus. The thick black solid line corresponds, in each 
plot, to a reheating temperature exceeding the thermal freeze-out temperature for all models along the curve, namely 
it gives the models matching the cosmological DM density we would also obtain in the standard picture without 
non-thermal DM sources: In the Mi-/i plane this happens, starting at small neutralino masses, close to the diagonal 
Ml = fj, since it requires a tuning of the right amount of Higgsino and Bino component in the LSP, suppressing the 
large Higgsino annihilation cross section with the Bino one, which in Split SUSY is extremely small. Would we have 
allowed for lighter sfermions and other light Higgs states, this curve would have moved only slightly further away for 
the diagonal, except when sfermion coannihilations or S-channel resonant annihilations on a Higgs take over in setting 
the effective thermally averaged annihilation cross section, as happens, e.g., in portions of the mSUGRA parameter 
space, see, e.g. [i^ - we will not discuss these exceptions here. As already mentioned, at about 1.1 TeV a pure 
Higgsino saturates the thermal relic density bound. Turning to the M2-IJ- plot, Winos have an even larger annihilation 
cross section than Higgsinos and the thermal relic density curve just goes from a pure Higgsino to a heavier pure 
Wino through a transient with large Higgsino- Wino mixing. Finally the behavior in the M1-M2 plane is more peculiar 
since from the structure of the neutralino mass matrix, Bino and Wino do not mix and, below the mass scale for 
a pure Wino thermal relic candidate, the tuning here is between the mass spitting between the Bino LSP and the 
second lightest neutralino and the lightest chargino, which are Wino-like and and whose coannihilations in the early 
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FIG. 4: Models with relic density equal to the central value from the 7-year WMAP data analysis, i.e. Q^h^ = 0.1123, in 
two-dimensional slices of the Mi, M2, parameter space, corresponding to the limit in which the third parameter is heavy (set 
to 10 TeV in the numerical computation), for a few values of Trh as indicated in the plots, and two values of Bx, namely 0.25 
for left panels and 10~^ for right panels. The filled areas correspond to regions violating the LEP lower bound on the chargino 
mass. 

Universe set the thermal relic abundance of a Bino LSP (there are chargino and neutralino coannihilations even for 
pure Winos and Higgsinos, but with less dramatic effects). Turning on the non-thermal component from the modulus 
decays, when Bx is large, essentially one just sees in the plots the scaling sketched in Eq. with Higgsinos and 

Winos saturating the WMAP preferred value for Vl^ with a progressively larger {av) as Trh decreases, and hence for 
a progressively smaller LSP mass, covering the whole parameter space for Higgsinos lighter than 1.1 TeV and Winos 
lighter than 2.4 TeV (as the Higgsino mass approaches the W boson threshold the cross section stops increasing; 
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this explains the shape of isolevel curves in that region). A detection at an accelerator of such LSP configurations, 
hopefully combined with a DM detection signal, would indeed be an indication of a non-standard cosmological phase 
at DM generation, with non-thermal production as primary scenario (there would also be further possibilities, such 
as, e.g., the increase of the Universe expansion rate at freeze-out induced by a quintessence component [43, |45| or a 
modification of the gravity theory (1^). When Bx is small, there is a smooth transition from the regions where the 
scaling in Eq. ([T8|) applies to those where annihilations stop playing a role and Eq. ([20l) applies instead; the latter 
makes even pure Binos, which, we underline again, in our sample MSSM setup have extremely small annihilation cross 
sections, become cosmologically viable, another configuration which, if singled out at accelerator and/or DM searches, 
would point to a non-standard early Universe cosmological history (in plots the filled region stands for the region in 
which the LEP bound on the chargino mass > 103.5 GeV is violated; Tevatron and the recent LHC constraints, 
such as [13,14^, are not shown since we have made just schematic assumptions on sfermions and not discussed at all 
gluinos, the particles most critical for a early discovery at a hadron collider). 

V. NON-THERMAL DM PRODUCTION IN THE G2-MSSM 

As an example of framework in which we can make definite predictions for the spectra of both the set of particles 
{Xa} E^nd the cosmological moduli, we discuss the case of the G2-MSSM Within a specific class of string 

M-theories, in this scenario the compactification of extradimensions gives rise to a TV = 1 supergravity in which SUSY 
breaking, due to the dynamics of a hidden sector, is trasmittcd to the visible sector by a combination of gravity 
(dominant contribution) and anomaly mediation. We briefly summarize here the main features of the spectrum, 
following Ref. [12]: In the G2-MSSM the visible sector can be described by a GUT theory broken into the MSSM 
at the unification scale Munif, at about 10^^ GeV, coinciding with the compactification scale. The RGEs boundary 
conditions are mainly functions of the gravitino mass TO3/2, which can be estimated from the UV theory parameters to 
lie in the range between ten and several hundred TeVs. The gauginos are expected to be the lightest SUSY particles; 
at Munif the gaugino masses are generated from a universal loop-suppressed gravity mediation contribution combined 
with a non-universal anomaly mediation term. The ratio of the gaugino masses to the gravitino mass depends almost 
linearly on the quantity 6 that parametrizes a threshold correction to the unified gauge coupling; in the following, S 
will be kept as a free parameter. The value of the masses at the electroweak scale is computed following the RGE 
evolution, including threshold corrections, such as the very large correction coming from higgs-higgisino loops, which 
is proportional to /i [49] ; whether the lightest neutralino is the Bino or the Wino depends on the sign and magnitude of 
this latter correction. For /i > and for 6 in the range — 10 < 5 < the lightest neutralino is a pure Wino, with mass 
in the range between about 100 and few hundred GeV (even the gluino is fairly light, mg < 1 TeV, a feature implying 
a rather rich phenomenology at LHC (sdHHI and making the model testable in the near future). Concerning the other 
states in the MSSM spectrum, the Higgsino mass parameter n and soft SUSY breaking term are generated by 
a Giudice-Masiero mechanism and are heavy, of the order of Sfermions are also heavy with a fiavor universal 

contribution to their soft masses at Munif being about TO3/2; RGEs affect mostly the third generation of squarks with 
the stops and the left handed sbottom becoming the lightest sfermions at the Electroweak scale (the left-handed stop 
mass becomes about 0.25 • m^/2^ the right handed stop and the left handed sbottom masses about 0.5 • 7713/2). In the 
Higgs sector tan /3 is fixed by /i and _B/i through the electroweak symmetry breaking condition and takes a value of 
order one; the light CP even Higgs is Standard Model like, while all the other Higgs bosons are heavy, again at about 
the r7i3/2 scale. In summary, the features relevant to discuss non-thermal DM production in this model are: The pure 
Wino LSP as DM candidate, as enforced by the proper choice of 5 (we will restrict to values 5 < —3 since they are 
theoretically favored [19]) and required to provide an annihilation cross section sufficiently large for the model to fit 
into the scenario in which the branching ratio for the decay of the moduli into the LSP is unsuppressed; A charged 
Wino as next to lightest SUSY particle, with a tiny mass splitting with respect to the LSP, about 200 MeV, due to 
the one- loop electroweak corrections to the neutralino and chargino masses [50| : The possibility for the moduli to 
decay into gluinos and third-generation squarks. While other SUSY particles do not play a role in our analysis, the 
relevant part of the spectrum will be computed here implementing the appropriate one-loop RGE running. 

For what regards the moduli fields, as already mentioned, in a string framework like the G2-MSSM, they arise in the 
effective supergravity theory after the compactification of the extra-dimensions. The theory predicts the presence of 
a large number of moduli fields with mass of order or heavier than the gravitino mass. In our numerical computation 
of the DM relic density, we follow the scenario outlined in Ref. ^25] : The modulus sector is composed by -|- 1 fields 
with N = 0(50 - 100) including: 

• 1 heavy modulus Xn with mx^ = 600 7713/2; 

• 1 meson field $ with 777$ — 1.967773/2; 
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FIG. 5: Left panel: Wino relic density versus the gravitino mass 7713/2 for a few sample values of Dx^ and S = —3; the upper 
horizontal scale shows the corresponding value of for for this specific value of 5. The band gives the 1 — a determination of 
the DM relic density from the 7-year WMAP dataset. Right panel: models with relic density equal to the central value from 
the WMAP data, in the plane versus 171^/2, obtained by varying 5 in the range — 10 < 5 < —3 and few sample values of 
Dxi ■ The filled area marks the region violating the LEP lower limit on the chargino mass; the region of the plane above the 
dashed line would correspond to models with 5 > —3. 



• N-1 light moduli Xi with mx^ = 1.96 7713/2- 

Their decay rates can be written in the form given in Eq. ([T]), i.e. it is proportional to the mass of the modulus 
to the third power and inversely proportional to the square of the reduced Planck mass. The constant in front can 
be computed esplicitly [1^ in the model; we keep -Dxjv for the heavy modulus and for the meson fixed to the 
benchmark values of, respectively, 2 and 710, while Dx^ , which is one the quantities the LSP relic density is mostly 
sensitive to, will be treated as a free parameter allowed to vary in the range between 4 and 16 (preferred range in the 
scenario considered here With this choice, the lifetimes for the three type of states is split to about: 10~^° s 

for Xn, 10^^ s for $ and 10~^ — 10~^ s for the light moduli. The branching ratios into SUSY particles of the decays 
are also calculable in this theory, with the main channel being into squark pairs, mostly the lightest stop, which in 
turn cascade down to the LSP and the chargino; in general, the branching ratio of decay of the light moduli into 
Susy particles is 25% with on average two DM particles produced at the end of the decay chain, giving Bxi ^ 0.5. 
Gravitinos are produced in the heavy modulus decay, while gravitino pair production in the decay of the meson 
and the light moduli is kinematically forbidden for the given values of and mxi , a choice quite natural for the 
G2-MSSM but still not totally general [ll|. We will comment further on this point below. Gravitinos are also long 
lived: 

1 "7i^/2 

producing one SUSY particle per decay, cascading again into one DM particle. 

The system in Eq. p4)) is solved numerically for this G2-MSSM setup, with the set of moduli just outlined and 
having chosen N = 99 and including the gravitino as V' field. The quantities which are kept as free parameters in our 
analysis are the gravitino mass, the parameter S (allowing to shift the ratio between LSP and gravitino mass ) and 
Dxi ■ Except for gravitinos, all other decay products are, for the moment, treated as particles in kinetic equilibrium. 
The system is evolved starting with the oscillations of the heavy modulus, when its initial energy density l/2m^^Mpi 
is equal to the radiation energy density, while the other moduli are included in the system at the beginning of their 
oscillations. For what regards the generation of DM, the relevant production phase is only the one from the decays 
of the light moduli, since the thermal DM component, as well as those from the decay of the heavy modulus and 
the meson, get diluted in the entropy injection phases. The dependence of DM conioving number density on 
temperature in this scenario is perfectly specular to those shown in Fig. [1] for models whose number density follows 
first a phase of the quasi-static equilibrium and then reannihilation. Gravitinos play a marginal role: produced in the 
heavy modulus decay, they get diluted and decay at a late stage (possibly after the end of the reannihilation phase 
for x) when Y3/2 is tipically 3 to 4 orders of magnitude smaller than the final y^, hence not contributing significantly 
to the DM relic density. 

In left panel of Fig. [5] we plot the neutralino relic density versus 77^3/2, for a few values of Dxi and a sample value 
for S, showing also on the upper horizontal scale the corresponding value of the Wino LSP mass for such given S. 
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Given that the reannihilation regime apphes, from Eq. (|18p we expect to be proportional to and inversely 

proportional to {av) and Trh, with the latter in turn approximately proportional to -D^ ^ ("^3/2)''^^5 see Eq. ([3]) where 
the scaling in the modulus masses has been replaced by the scaling in terms of the gravitino mass. In the limit in 
which the Wino pair annihilation cross section just scales with the inverse of the square of the Wino mass, one would 
find: 



^2 [F{5)f (m3/2 



3/2 



^xh' (X ' (27) 

where the function F{5) parametrizes the quasi-linear relation between and m3/2- In the plot, the result of the full 
numerical solution roughly confirms these approximate scalings, except for small for which {av) is not inversely 
proportional to m^. To match the experimental value the DM abundance, lighter TO3/2 and larger Dxi are favored. 
In the right panel of Fig. [5] we consider the plane versus and, varying 5 and for a few values of , we plot 
models that have ri^^^ equal to the mean value from the WMAP data; the plot illustrates the fact that, even in a 
model as constrained as the G2-MSSM, there is still a rather large sensitivity to the parameters setting the theory 
at high energy. A relic density compatible with cosmological measurements is obtained for LSP lighter than about 
300 GeV and for reheating temperatures in the range between about 100 MeV and 1 GeV. The results of our analysis 
are consistent, as an overall picture, with the results presented in Refs. (25l. [52j. although there are slight numerical 
differences when comparing model by model; most likely these differences stem mainly from the determination of the 
mass spectrum of the G2-MSSM which is probably less accurate in our work, although the more careful numerical 
treatment implemented here for the relic density calculation may have some impact as well. As a final remark, we 
mention that we have also crosschecked the result that, to obtain a relic density compatible with the DM density 
as measured by WMAP, it is necessary to forbid the decay of the light moduli into gravitinos; in case it is not, in 
all G2-MSSM setups, gravitino decays become the main dark matter source, at a stage when reannihilations are 
inefficient, largely overproducing dark matter. 



VI. KINETIC EQUILIBRIUM AND DECOUPLING IN THE G2-MSSM 

In Section [n] we have emphasized that Eq. ([7]), tracing the evolution of the number density of the Xa particles, has 
been written assuming that kinetic equilibrium between the Xa states and the thermal bath particles is maintained at 
all stages over which the comoving number density changes. Whether this assumption is valid or not depends on the 
efficiency of the scattering processes on thermal bath particles, an issue with is usually addressed invoking crossing 
symmetry arguments relating the scattering to the annihilation cross section; in most explicit models however the 
two processes are not related via crossing symmetry and one should actually study this problem case by case. We 
focus here on the G2-MSSM (a more general framework with non-thermal Wino DM will be also considered at the 
end) and discuss the steps which should be followed when relaxing the hypothesis of kinetic equilibrium, introducing 
a more general set of Boltzmann equations. 

The energy spectrum of SUSY particles produced in the decay of moduli is usually very different from the thermal 
distribution; in particular in the G2-MSSM scenario, light particles are generated in the decay of very heavy fields. 
The cascade process generally starts with the production a pair of squarks, followed by their decay into gluino and 
quark, and with the gluino in turn decaying with a three body process into the LSP, the Wino-like chargino or the 
Bino, with branching ratios depending on parameters in the model. As a last step, the Bino decays as well into the 
chargino or the LSP, while the chargino, given the small mass splitting with respect to the LSP, has a longer lifetime. 
The chargino decay occurs either through a two body process in which a pion is produced together with the LSP, or 
through a three-body in which a neutrino and an electron are produced; the rates of these processes are given by, 
respectively, [53.. ,54]: 

rx±,2b = -^Am2y^Am2 -m2 and r^± ,3^ = , (28) 

where Am^ is the chargino-neutralino mass splitting, /^r — 93 MeV is the pion decay constant and Gp is the Fermi 
constant. The two-body decay is dominant when kinematically allowed; this is the case in the G2-MSSM, since the 
minimum mass splitting between charged and neutral Wino, induced but electroweak radiative corrections to the two 
masses is Am^ ~ 160 MeV [HO] • We have studied the decay chain of the moduli with the package PYTHIA [55| for a 
few sample benchmark models in the G2-MSSM, assuming a stable Wino-like chargino, and found energy distributions 
for the Wino-like neutralinos and charginos which are typically peaked at E/m^ ^ 10 and with very broad tails up 
to the kinematical threshold; among decay products, the number of charginos is typically about 3 times larger than 
the number of neutralinos. 
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The injected ultra-relativistic particles lose energy via scattering on thermal bath states. Were these processes 
inefficient, the non-thermal DM generation would give rise to a model of the Universe with warm or even hot DM, a 
possibility which has been investigated, e.g., in Refs. [5B,[53|. As a first rule of thumb, the energy depletion is efficient 
whenever the relative energy loss rate times the time interval the over which the process is active, which we indicate 
as At, is larger than 1: 

1 dE\ , , , 

In our case this condition needs to hold from the relativistic regime down to the non-relativistic low-temperature 
environments induced by the reheating phase. The expression for —dE/dt is in the form: 



(30) 



where r(i?) the scattering rate for the process under scrutiny, integrated over the phase space distribution functions 
of the thermal bath particles in the initial state and the phase space of the out-scattered particles. The expressions 
we will report below are derived in the limit of small momentum transfer between the non-thermally produced states 
and the thermal bath particles; the latter on average have energies equal to about 3 T. The small momentum transfer 
approximation holds whenever the non-thermal particles are non-relativistic in the CM frame of the scattering process, 
namely for [56j : 



m\>&TE. (31) 



Assuming instantaneous production at reheating, this relation can be translated into: 

10 



■ lOOGeV/ \E/m 

a condition which is satisfied in the region of the G2-MSSM parameter space providing a viable DM candidate. 
Charginos lose energy via electromagnetic interactions and their energy loss rate takes the form [58, 59]: 



dE\ wa^T^ . 

— ^A (33) 



- 



dt J 



with A ~ 0(1). The elastic scattering of a Wino-like neutralino on a background lepton is very inefficient, since 
it proceeds via a Z boson or a slepton exchange and the corresponding amplitudes are suppressed, respectively, by 
the tiny higgsino fraction in the LSP and by the slepton masses, which in the G2-MSSM are very heavy. Whenever 
kinematically allowed, the dominant effect is the inelastic scattering into the charged Wino, which is mediated by a 
W boson. There are then two effects making a neutralino produced in the decay of the modulus lose energy, namely 
the energy loss in the inelastic scattering itself and the fact that the produced chargino will efficiently lose energy. For 
relativistic neutralinos, the inelastic scattering rate and the energy loss rate in inelastic scatterings are, respectively, 
given by: 

(a,b) 



X \ "-x 



dE\ 's—^ IQgwabGp f ■m^Am^\ 5 / E \ ^ f o^"^ 



dt / ,.o TT"^ V 2ET J \m^j \ m 

^ (a,fe) 

where, considering the generic process in which the heat bath particle a is scattered into the particle h via an interaction 
vertex with a W boson, we have included in the coefficient gwab the product of the number of internal degrees for a, 
that for 6, as well as a rescaling factor in case the coupling constant in the vertex is different from the SU(2) gauge 
coupling g (e.g., for the scattering process + ^ + ^e, gwab = 8); the sum goes over any (a, h) thermal bath 
particle pairs. 

In Fig. [S] we consider two of the G2-MSSM singled out in the previous Section as models embedding a viable DM 
candidate, at the light and heavy ends of the mass range displayed in Fig. O i.e. two models with Wino masses, 
respectively, of 103.5 and 300 GeV, obtained for = 16, (5 = —3.5 and 6 = —3 and gravitino masses of 107 and 
460 TeV, and corresponding to scenarios with approximate reheating temperatures of 100 MeV and 900 MeV. For such 
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(a) Case \:my^ = 103.5 GeV, Tr = 100 MoV, E/m^ = 1.005. (b) Case 2:mx = 103.5 GeV, Tr = 100 MeV, E/m^ = 10. 
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(c) Case 3;m^ = 300 GeV, Tr = 900 MeV, Elm^ = 1.005. (d) Case i-.m^ = 300 GeV, Tr = 900 MeV, E/m^ 
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FIG. 6: Ratios between the effective rate of energy loss rate —1/E ■ dE/dt (dashed lines), or of the scattering/decay rate F 
(solid lines), to the Universe expansion rate H, for a few processes involving charginos and neutralinos. The upper panels refer 
to a G2-MSSM DM model with = 103.5 GeV and Trh = 100 MeV, while the lower panels to one with = 300 GeV and 
Trh = 900 MeV; the plots on the left hand-side refer to non-relativistic particles, E/m^ = 1.005, while those on the right-hand 



side correspond to a sample relativistic case, E/m^ 
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models we plot ratios of scattering and decay rates F, or of relative energy loss rates —1/E ■ dE/dt, to the Universe 
expansion rate H; in the panels on the right-hand side, results are shown for relativistic particles, E/m^ = 10, while 
on the left-hand side the non-relativistic limit is considered, E/m^ = 1.005. To sketch the efficiency of the chargino 
energy losses, the appropriate timescale Ar in Eq. is the shortest between the chargino lifetime and the timescale 
for back-scattering of the chargino into the neutralino, i.e. the rule of thumb condition in Eq. (|29p holds whenever 
the curves in plots corresponding to the chargino energy loss lie above the curves for the decay rate and the inelastic 
scattering rate. More quantitatively, for the two processes, these ratios are: 



1.86- IQ-^A 1- 



'1/2 



T 



IMeV 



X , 1— scat. 



6.77 • 10^2^ (^^^j ( 



100 GeV/ \E/m 



^ ^leOMeVy 1^100 GeV 
Amv / 
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(36) 
(37) 



with the smallest of these being much larger than one in all cases except for relativistic charginos injected at tem- 
peratures of the order of 1 GeV or larger. The latter is however the regime in which inelastic scatterings turning 
a neutralino into a chargino and viceversa are extremely efficient (as shown in the plots the rate for this process is 
many orders of magnitude larger than H) and the energy loss rate via this process is also very large (the relevant 
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timescale is now ^ H~^). This shows that the relativistic charginos injected at any of the temperatures of interest 
in our model instantaneously thermalize. For what regards neutralinos, in the relativistic limit, the energy deple- 
tion is guaranteed by inelastic scatterings and by chargino energy losses down to background temperatures of about 
2 MeV; however when becoming non- relativistic and at low temperatures, the rate for inelastic scatterings becomes 
smaller than H and the assumption of kinetic equilibrium may not hold any more. To study the evolution of the 
system at low temperature and model kinetic decoupling, we follow the approach of Bringmann and Hofmann [l3l | 
(see also 60]) who have developed a formalism to treat kinetic decoupling starting from the Boltzmann equation for 
the phase-space distribution function of the WIMP DM candidate; we extend here their treatment to the case of two 
co-annihilating particles. Let fy.o(jp,t) and fj^±{p,t) be the phase space distribution functions for, respectively, the 
neutral and charged Winos. We have just shown that charginos are kept into kinetic equilibrium at all temperatures of 
interest for your problem, so we can assume that the shape of the chargino distribution traces the thermal distribution 
function, namely: 



/x± (p, t) cx fll {p, t) 



(38) 



On the other hand, the distribution function of the neutralinos could have a shape which is slightly different from the 
thermal one, since we have shown that energy losses may not be efficient in the non-relativistic regime; this departure 
is parametrized defining the temperature of neutralinos T^o through the second moment of the distribution function: 



(2^ 



9x" fx° (P' = T^o (t) n^o (t) . 



(39) 



For neutralinos in kinetic equilibrium, T-^o coincides with the thermal bath temperature; after kinetic decoupling the 
neutralino temperature will scale instead as T^o oc . 

The two distribution functions obey the system of coupled Boltzmann equations: 



{dt ~ Hp-Wp) f^a{p,t) = lc^o[/^o,/^±] 



(40) 



idt-Hp-Vp)f^±{p,t) = -C^±[/^o,/^±], 

where C stands for the collisional operator, embedding all interactions involving neutralinos and charginos, namely 
annihilation and scattering processes, as well as the production of neutralinos and chargino from moduli decays and 
the neutralino source from chargino decays. Integrating these equation over phase space one obtains two equations 
for the time evolution of the neutralino and chargino number densities: 



dt ^ 



dn^± 

^ '3Hn^± 



dt 



g^on^± — g^±n^o exp 



eq eq 



F^o^^i + F^± 



g^±n^o exp 



At 



T 



eq eq \ . \ ^ ^ 



Arn^ 
T 



(n-o)^ 



(41) 



In these equations, the first term on the right hand sides accounts for inelastic scatterings of neutralinos into charginos 
and decays of charginos into neutralinos (as well as the inverse processes); g^o and g^± are the number of internal 
degrees of freedom for the neutralino and chargino, while F^± is the chargino decay rate as obtained including the 
two contributions in Eq. (|28| . For inelastic scatterings we have assumed that: i) the diagram with W boson exchange 
in the t-channel is the dominant one (since those with sfermion exchanged are highly suppressed) ; ii) the momentum 
transfer in the t-channel is small and the collision term can be computed expanding in its powers, see also fl^, IgQ; 
in) Am^ and T are small with respect to and only the lowest order terms in an expansion in Am^ and T give 
sizable contributions; under these assumptions, we find: 



y-' gWab 
(a,f,) 



■T3 (Ar 



12T2) 



(42) 



(some further details and a sketch of the derivation of this expression is given in Appendix 1^). When including pair 
annihilation terms in Eq. (I4ip we have taken advantage of the fact that it involves non relativistic particles annihilating 
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mainly via S-wave processes and hence the cross section has a very mild dependence on momentum, allowing then us 
to write an expression analogue to thermal case also when the neutralino distribution function starts deviating from 
the kinetic equilibrium value. Furthermore, since the relativistic particles injected from the moduli decays mostly 
lose energy as charginos, decaying afterwards into neutralinos, we have simplified the treatment including these as 
a source function of "thermal" charginos only. Obviously, summing the two equations one retrieves Eq. ([7]) with 
being the sum of the number density for the two coannihilating species. 

Taking the second moment of the first equation in the system in Eq. (j40p , one finds that the neutralino temperature 
T^o obeys the relation: 



(the derivation of this equation is also sketched in the appendix) . 

The numerical solution of the problem proceeds now analogously to what done so far. After the appropriate change 
of variables, the system in Eq. (|4T|) replaces Eq. ([7]) in the system of Eq. ([T4|) . The explicit solution for n^o[t) and 
n^± (t) are then implemented in Eq. (|43p to find T-^o (t) . 

Our first application is to the G2-MSSM models singled out in the previous Section as cosmologically favored. As 
we had guessed in the analysis we performed at the level of energy loss and scattering rates and shown graphically in 
Fig. |6l the departure from kinetic equilibrium tends to be at a temperature sensibly lower than the nominal reheating 
temperatures for these models (which are of the order of 100 MeV or larger). The numerical solution indeed shows 
that the ratio n^± /n-^o tends to follow very closely the ratio of the thermal equilibrium number densities / n^o over 
the whole phase of DM production in the moduli decays, as well as at later times. The solution of the equation for 
the neutralino temperature shows that kinetic equilibrium is maintained up to a temperature of the order of 10 MeV, 
independently of the neutralino mass since, in the non-relativistic limit, the inelastic scattering rate (which together 
with chargino electromagnetic interactions enforces the equilibrium) depends only on the chargino-neutralino mass 
splitting which is essentially the same over the whole range of selected models. The transition between T^^^ = T to 
the regime T^o oc takes place on relatively short timescales; since at 10 MeV non-thermal production has become 
irrelevant, we would have found the very same scaling when computing the kinetic decoupling for a population of 
thermal particles: the evolution of the number density ratio and of T^o for G2-MSSM DM models are those shown as 
black lines in Fig. [7] and labelled, respectively, 'thermal distributions' and 'standard decoupling'. 

To illustrate the impact of non-thermal production and non-standard cosmologies on the kinetic decoupling process, 
we allow then for a slight variant to the underlying particle physics framework, still referring to a pure Wino as DM 
candidate but assuming now that the reheating temperature can be reduced to values much closer to the bound from 
BBN than in the G2-MSSM. In Fig.[7]one sees a modification with respect to the standard case when the gap between 
reheating temperature and standard kinetic decoupling temperature is reduced, i.e. for Trh equal to about 20 MeV or 
lower: The additional DM source makes the ratio n^± /n-^o differ from the ratio of thermal distributions. The impact 
on T-^Q is two folded: the chargino decays tend to populate the system with neutralinos that are on average more 
energetic than for a thermal distribution, delaying the onset of the regime T^^o oc and making the transition into 
this regime to be less sharp; at the same time, if Trh is so low that reheating increases significantly the expansion 
rate of the Universe H at the time of kinetic decoupling (Trh = 8 and 5 MeV in the plot) the departure from Ty,g — T 
tends to be anticipated. This latter feature was already pointed in [61i] . showing that the non-thermal production 
could induce higher kinetic decoupling temperatures compared to the standard case; in case of Wino DM, however, 
the production and decay or charginos in the moduli decay has always a larger impact. The kinetic decoupling 
temperature is directly related to the minimum mass scale for structures in the Universe; we have shown here that 
even in case of injection of particles with efficient energy losses on the thermal bath, the low-temperature non-thermal 
production can leave an imprint on structure formation. The development of a precise numerical treatment of the 
kinetic decoupling is then a valuable tool to test this class of models. 

Finally, in the examples considered here, we find a marginal change in the DM relic density when computing it in 
the case when we trace the the number densities of the individual coannihilating species as opposed to the case when a 
single equation for the sum of number densities is solved; this is due to the fact that the departure of the ratio n^± / n-^o 
from the ratio of thermal distributions takes place only when such quantity is very small (moreover, in our examples, 
the annihilation rates for each of the coannihilation channels are comparable) . Considering however models for which 
crossing symmetry arguments between annihilation and scattering cross sections are even more severely violated, one 
should find cases in which the standard thermal assumption is invalid at higher temperatures, possibly even close to 
the chemical freeze out temperature; in those cases there should be a sizable change in the relic abundance as well 
and the formalism we developed would be suitable for an accurate computation of the relic density for such case. 
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FIG. 7: Left panel: ratio of the chargino number density over the neutrahno number density for several values of Trh- Right 
panel: Ratio T^„/T as function of the temperature of the thermal bath for same values of Trh- Plots are obtained for a Wino 
with mass equal to 200 GeV, however results depend only on the chargino-neutralino mass splitting which is about 160 MeV 
in the scenario under consideration. 



VII. CONCLUSIONS 



Non-thermal dark matter production is a viable alternative to the standard paradigm of WIMPs as thermal relics. 
It is a well-motivated scenario arising naturally in several particle physics frameworks, including SUSY standard model 
extensions within supergravity and superstring theories. Moreover, an epoch of entropy injection at an intermediate 
phase between the reheating at the end of inflation and the onset of BBN dilutes dangerous long-lived relics, such 
as gravitinos. We have reviewed how to introduce a system of differential equations to treat a generic case of non- 
thermal dark matter generation and implemented an efficient and accurate numerical scheme for the computation; 
such scheme has been interfaced to DarkSUSY numerical package and will be released together with an upcoming 
version of the code. The interest in this scenario has been recently boosted by the fact that, generically, it favors 
dark matter candidates with pair annihilation cross sections larger than in the thermal WIMP framework, possibly 
suggesting a connection with the very large annihilation rates which would be needed to explain with a dark matter 
induced component recently detected anomalies in cosmic-ray lepton fluxes, such as the rise in the positron fraction 
measured by the PAMELA detector; we have considered, within a toy model, what range of reheating temperatures 
would follow from such identification. Focussing on SUSY models, we have discussed the impact of this non-standard 
cosmology in selecting the preferred mass scale for the lightest SUSY particle as dark matter candidate, an issue with 
a direct impact, e.g., on the interpretation of new physics eventually discovered at accelerators including the LHC. 
In the second part of the paper we have concentrated on a more predictive model, the G2-MSSM, and questioned in 
further details the underlying assumptions in the standard solution of the Boltzmann equation for the dark matter 
component; in particular, we discussed how to verify whether kinetic equilibrium holds along the whole phase of dark 
matter generation, as well as the validity of the factorization usually implemented to rewrite a system of coupled 
Boltzmann equation for each number density of a set coannihilating particles as a single equation for the sum of all 
the number densities. As a byproduct we developed here a formalism to compute the kinetic decoupling temperature 
for a system of coannihilating species, which can be applied also to other particle physics frameworks, also in case 
standard thermal relics within a standard cosmology. 
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Appendix A: Evaluation of the collisional operators 



In this appendix we will sketeh how to compute the operators C as introduced in the system of coupled Boltzmann 
equations (1401) . As a sample term we discuss how to deal with the contribution to the neutralino collisional operator 
coming from inelastic processes of the type: 

x''{P) + a{K)^X^{P') + b{K') (Al) 

where a and b are thermal background particles, with four-momenta, respectively, K = (fc,k) and K' = (fc',k'), 
while P = {E, p) and P' = {E' , p') denote the four- momenta of neutralino and chargino. Summing over all avalaible 
thermal bath pairs (a, 6), such contribution to the collisional operator, normalized to the neutralino energy E, takes 
the form: 

E J (2iT?2kJ (2iT?2k' J (2TTf2E' 2E ^ ' ^ ' 



(a.b) 



■ [h{k'){l ~ famh^ip') - faim - Mk'))f^oip)] ; (A2) 

where is the modulus squared of the scattering amplitude, averaged over the initial spin states and summed over 

the final spin states. In the following we assume that a and b are massless and described by Fermi-Dirac distribution 
functions fa{k,t) and fb{k',t) (to shorten the notation the indices a and b will be dropped). The exchanged four- 
momentum is indicated as (w, q) = (E — E' , p — p'); for kinematical reason, the transferred momentum is constrained 
to be of the order of the heat bath temperature and is small compared to the initial energy and mass of the neutralino. 
The ratio between the mass splitting chargino-neutralino and the neutralino mass will be also assumed as a small 
parameter. Furthermore we have initial conditions such that neutralinos have nearly thermal distributions, implying 
V ^ y^T/mx ^ 1, being the neutralino velocity defined as v = p/ E. Under such conditions, we can just take the 
non-relativistic limit of the collision term and, in addition, expand it respect to the quantities T/m^ and Am^/m-^. 
Using these informations we can eliminate the dipendece on p' by Taylor expanding f~^±{p') as: 



4± (p') c (p) - q . Vp4± (p) + i(q • Vp)^/^± (p) + ... (A3) 



This allows us to freely integrate over the three-momentum component of the delta function in Eq. (jA2[) . Using now 
the relation: 

/(fc')(l - f{k)) = exp (-^) /(fc)(l - /(fc')) ^ exp (-^) /(fc) (A4) 

with: 

e^J^^)c.e.p(^) fl-^!!^-q:Z + ^+(q-")(q-")) (AS) 



we need to compute: 

Cx'>,is[fx°] _ ^ gwabQx^ f d^k ^^^^ f d^k' 

{a.b) 



E 



y ~9wab9x^ f f d^\M\l, 5 {E' + k' - E - k) 

^ 2b&T:'^ EE' J k ■'^ ' J k' ^ ^"'^ ^ ' 



q ■ V dfx^ , q V 

. y 2T • V dp T -""^ 

l(q^ 3(q.v)(q.v)V/x±g^ 



2 \ u J dp 



(A6) 



which is the analogous to the expression for the expansion of the collisional operator obtained in [13| . As an example 
we sketch the calculus of the first term in the square bracket. The invariant amplitude takes the form: 

\M\l, = 64G| {{PK){P'K') + {PK'){P'K) - (m^ + Am^)m^KK') (A7) 
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At the leading order in T /m^ and ^m^jm^ we can write: 



(^2{k - V k)(fc' - V k') + ^(a;2 - q^){l - v^)j + y/l-v^ J^- 



, , Amy , , , Aruy 
-fcv • k' ^ - k'v ■ k ^ 



, Ar 



^,k^^k^^2v.kv.k'^ + 

3„2^(v•q) 



(A8) 



The first step of the integration is an average over the directions of v. It can be done by use of the foUowing results: 

f^SiW-^r■^) = —9{v^q^^W^) (A9) 
J 47r 2vq 



6{W - 



2vq q^ 



f ^S{W - V • qyv^ - T^O{v^q^ - W- 
J 47r 2vq 

f ^5(W - V • qyv^v'' = —9(v^q^ - W^) 
J 47r 2vq 



2g2 



2g4 



- Sv^q^W 
2^ 



-q'-q^q'' 



where we defined W = u + \/l — v'^Am^. The further 3 integrals in Eq. (IA8|) are most easily performed in in the 
variables uj, q and k. Actually the last step, the integral in the variable k needs to be performed numerically; an 
analitic expression, which traces rather accurately the numerical result, can be obtained by replacing the Fermi- 
Dirac distribution f{k) with the exponential scaling exp(— /c/T). At the first order in T/m^ and Ato^/to^, also 
remembering that, under our assumptions, Irr? ~ T jm^ this gives: 



^/(fc) / '^^\M\l,KE'^k'-E-k) = 



(Ato^ + eAm^r + 12T2) ( 1 



-2AmlT^ (Amy + 2T) — + -T (Amt + SAto^T + 32ArnlT^ + 114ATOyr3 + 144T*) — , 



+ 

2 1 



e— .(AlO) 



After computing all integrals, the inelastic scattering contributions to the neutralino and chargino collision terms are 
then found to be: 



E 



E 



2Gpgwabgxi 



AT^ (Ami + QAfUyT + 12T^) 1 



Am^ 



2Am^T2 (Am^ + 2T) — + -T {Am\ + BAm^T + i2Am\T'^ + 114AmxT3 + 144r^) ^ 



m^- 3 



/x± -/x«e" 



Am^T^ {Ami + ^^m^T + 12T'^ 
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P ■ Vp/x± 



- -T^ (Ato^ + lOAmi^T + GOAm^T^ + 240Amxr3 + 480r'*) ( Ap/^± 

-2r2Am^ (Am^ + GAm^T + 12T2) j 
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P ■ Vp/y± 



(All) 
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-2Am2r2 (Am^ + 2r) — + -T (Ato^ - 4Ato^T - lOAm^T^ + 30ATOyr3 + 144T4) Ar] f /^oe"^ - /^i) 



:Am,T3 (Am^ + 6Arn,r + 12T^) (^/,o + ^1^° 



^T''' (Am^ + lOAm^T + GOAm^T^ + 240Amxr3 + 480T'') ( Ap/^o + ^ ' ^^f^" + ^t-J^o ) 



-2T^Am^ {Ami + 6Am^T + 12r2) -^/xoe"^ }> (A12) 



m2 



The ordering of the terms is such that, when integrated over the momentum p of the neutrahno (first equation) or 
the chargino (second equation), the terms on the third and forth row of the two expressions cancel out while, for 
what regards the others, it can be seen that once summing the two equations one obtains a term proportional to 

— 3^ (^fxo^ ~ /x±) which cancels out too. The Boltzmann equation in (H^ are then obtained after the 

momentum integration p and using the fact that T/m^, Am^/m^ ^ 1 which allows to keep just the 0*^ order terms 
in both (IXTT|) and (|XT2|. 

The equation for the ncutralino temperature can be obtained from the second moment of the Boltzmann equation 
in phase space: 

/ (2^ ^"""^^ -Hp Vp) f^o ija) = J .9x« ^C^o [f^o] (A13) 

The left hand side can be rewritten as: 



(ITm dn^a 
3n^o —P- + 15HT^.n^o + 3T^o — ^ (A14) 
^ dt ^° ^ ^ dt ^ ' 



For what regards the right-hand side, we have that the contribution from annihilations, can be computed using the 
S-wave approximation. In this case in fact we can assume that the dipendence of the pair annihilation cross section 
on the momentum can be neglected and use the same factorization implemented when assuming kinetic equilibrium, 
i.e. : 

^s3 9l92 {(JV) (/1/2 - /f'/s'') = (ctw) (ni7l2 - ni^eqn2,eq) (Al5) 



(2^)^ (2^)^ 



with 



Jd^pid^P2i(Tv)flf2 

taking second moments of the distribution functions, one gets a the term: 

3M^oT^„ (^(crt))^o^o (n^o - «x°,oq) + i'^^)x°x^ {n^on^± - n^o^^qn^i^^q)) (A17) 
which cancels out against the term on the left hand-side proportional to dn^^/dt. 
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